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We develop a powerful numerical algorithm for calculat- 
ing the density of states p(ui) of the fluctuating gap model, 
which describes the low-energy physics of disordered Peierls 
and spin-Peierls chains. We obtain p(u) with unprecedented 
accuracy from the solution of a simple initial value problem 
for a single Riccati equation. Generating Gaussian disorder 
with large correlation length f; by means of a simple Markov 
process, we present a quantitative study of the behavior of 
p(ui) in the pseudogap regime. In particular, we show that in 
the commensurate case and in the absence of forward scat- 
tering the pseudogap is overshadowed by a Dyson singularity 
below a certain energy scale ui* , which we explicitly calculate 
as a function of £. 

PACS numbers: 71.23.-k, 02.50.Ey, 71.10.Pm 

The fluctuating gap model (FGM) describes the low- 
energy physics of one-dimensional fermions subject to 
static disorder potentials. The first quantized Hamilto- 
nian of the FGM can be written as |l| 

H = -iv F d x cr 3 + V{x)<7 + A(x)a + + A*{x)a- , (1) 

where V(x) and A (a;) are random potentials describ- 
ing forward and backward scattering, vf is the Fermi 
velocity (henceforth we set vf = 1), <J% are the usual 
Pauli matrices, fn is the 2x2 unit matrix, and a± = 
\{cr\ ± io-i)- Eq.(|l|) emerges as the effective low-energy 
Hamiltonian in different physical contexts. For example, 
fluctuation effects close to the Peierls transition in quasi- 
one-dimensional charge-density wave systems can be de- 
scribed by this Hamiltonian. In this case A[x) describes 
the time- independent part of the fluctuating Peierls order 
parameter, the probability distribution of which can be 
obtained from a Ginzburg-Landau expansion of the free 
energy H. For commensurate chains A(x) can be chosen 
to be real, whereas it is complex in the incommensurate 
case Truncating the Ginzburg-Landau expansion 

at the second order, A(x) is approximated by a Gaus- 
sian random process, with finite average (A(x)) = A av 
below the Peierls transition and A av = in the disor- 
dered phase. For commensurate chains the correlator 
of A(x) = A(x) - A av is (A(x)A(x')) = /^e-W-'Ve, 
whereas in the incommensurate case (A(x) A* (x')) — 
A 2 e -|*-x'|/« and (A(x)A(x')) = 0. Here £ is the or- 
der parameter correlation length, which diverges at the 
Peierls transition. The Hamiltonian ([[]) describes also the 
low-energy physics of disordered spin chains , which 
can be mapped onto disordered fermions by means of the 
usual Jordan- Wigner transformation. In many cases the 



filling of the effective fermionic system is then commen- 
surate with the lattice, so that A(x) is real. 

The fundamental quantity which determines the ther- 
modynamics of the model ([!]) is the density of states 
(DOS) p(lo). In general, one has to rely on approxima- 
tions to calculate p{u>) or its disorder average (p(tu)), but 
in special limits exact results are available. Besides the 
trivial case where V(x) and A(x) are constant, the exact 
(p(co)} can be obtained by various methods fs]-[7| in the 
white noise limit £ — ► 0, A s — > oo, with A^£ — * const. 
For real A(x) and V(x) = the average DOS is known 
to exhibit, for sufficiently small A av , a Dyson singularity 
H at u> = 0. In Ref. ||] we have shown that this singular- 
ity survives for arbitrary £ < oo. A recursive algorithm 
due to Sadovskii |l(| does not reproduce the Dyson sin- 
gularity, so that this algorithm cannot be exact. In fact, 
a subtle flaw in this algorithm has recently been found 
by Tchcrnyshyov 0. Because Sadovskii's algorithm (or 
generalizations of it) has been used in different contexts, 
e.g. to explain the pseudogap phenomenon in the normal 
state of the cuprate superconductors |]ll] , |l2l , it is also im- 
portant to investigate its validity for complex A(x) with 
quasi-long-range correlations. 

In this work we develop an accurate algorithm which 
allows us to investigate the regime A s £ ^ 1 where no ex- 
act solution is available. We find that in the commen- 
surate case the pseudogap is overshadowed by a Dyson 
singularity below a cross-over energy lo* , which we deter- 
mine as a function of the correlation length £. We also 
consider the incommensurate case for which Sadovskii's 
solution turns out to be qualitatively correct but leads to 
a wrong ^-dependence of the depth of the pseudogap. 

Riccati equation - In the following we use the special 
symmetries of the continuum model ([j]) to show that the 
DOS can be obtained without ever calculating the eigen- 
values of H JI3]. Instead, we obtain the DOS from the 
solution of a simple initial value problem for a Riccati 
equation. This will enable us to calculate p(co) with un- 
precedented numerical accuracy. For a given realization 
of the disorder the local DOS of the Hamiltonian (Q) can 
be defined by 10 

p(x,Lu) = -n^ 1 lm r Ti[a3G R (x,x,uj)] , (2) 

where the retarded 2x2 matrix Green function 
G R (x,x' ,u>) satisfies 

[id x - M(x, oj + iO + )\g R {x, x', ui) = 5{x - x')a , (3) 
M(x,lu) = [V(x) + A(x)a+ + A*(x)a-]o- 3 . (4) 

We now make the non-Abelian Schwinger-ansatz Pjlq] 
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g R (x,x' » = U{x,u)Gq{x - x')U- l (x' » , (5) 

where U(x,oj) is an invertible 2x2 matrix and Gq~(x) is 
the Green function to the operator id x + i0 + <73, i.e. 



'(x) o 
-6{-x) 



(6) 



In the following the w-dependence is suppressed. The 
ansatz (||) indeed solves Eq. (||) if [/ (x) satisfies 

[id x - M(x)\ U(x) = (7) 

with the boundary conditions 

£/ 12 (-oc) = C/ 2 i(oc) = 0. (8) 

Two different solutions of Eq. (0) are given by 

U+(x) = Texp[-z f M(y)dy] , (9) 

J — OO 

U-(x)=T- 1 exp[i M(y)dy], (10) 

J X 

where Texp is the path-ordered and T _1 cxp is the 
anti-path-ordered exponential function. Because M> — 
a^Mas and TrAf = 0, the matrices U a satisfy = 
(JsU^cts and AetU a = 1, which means that they be- 
long to the non-compact group SU(1, 1). Thus, the ele- 
ments of the U a satisfy U a22 — U* n , U a \ 2 = U* 21 , and 
If/ail | 2 — |J7 a 2i| 2 = 1. While each U a (x) only obeys one 
of the two conditions (ph. the combination 



U(x) 



_ J_ f E/_n(a:) U +12 (x) 

U-2l(x) U +22 (x) I ' 



(11) 



satisfies both boundary conditions. Here u = 
t/_n(— oo) = U +22 (oo), so that detU(x) = 1. Denot- 
ing the first column of the matrix U a by u a and the 
second column by v Q (so that v a = aiu*), we obtain 
from Eqs.@ and (|lj) 

rRt i \ fat ft i-WulW 
y [x, x , uj) = — i \ 0(x — x ) 

I 

v + (s)vt(s') \ 
+ 0(x - xj > (73 . (12) 



Here constitutes adjungation of u + , so that u is 
a 2 x 2-matrix. Equivalent but more complicated forms 
of Eq.(12) were first derived by Abrikosov and Ryzhkin 
ge§. Combining Eqs.(§) and @, we get 



1 U-uUln + U-nUjn 

P(x.uj) — —Re ■ - — 

' y ' it U-iiUXn-U-nUln 



(13) 



Since this expression only depends on the ratios <E> a (x) 
— iU* 21 (x) /U* 11 (x), we may also write 



1 1 + (*) 

P(X,LU) = — tie- . . ; - , . 



(14) 



Using Eq.(Q) we find that the $ a (x) are both solutions 
of the same Riccati equation, 

d x $ a (x) = 2iw{x)$ a {x) + A(x) - A*(x)<P 2 Q (x) , (15) 

where we have introduced uj{x) — lu — V(x). Similar Ric- 
cati equations have recently been obtained by Schopohl 
]l7| from the Eilenberger equations of superconductivity. 
To specify the initial conditions, let us assume that out- 
side the interval [0, L\ the potentials V(x) and A (a;) are 
real constants, Voo and Aoo- From the definition of $ a 
we find that Eq.([l^) should then be integrated with the 
initial conditions 



* +( o) = MD = Ji-^ + ^ ) (i6) 



where the square root has to be taken such that for 
Aqo — > the right-hand side of Eq.(|l6|) vanishes. Note 
that the initial values are simply given by the stable 
stationary solution of the Riccati equation (|l5|) with 
V{x) = Voo and A(x) =A W . 

The case of a discrete spectrum - For (uj — Voc) 2 < A 2 ^ 
the spectrum turns out to be discrete |18|: Introducing 
ifa(x) via $ Q (x) = e lVa ^ the phases satisfy 

d x y a (x) = 2u(x)-2\A(x)\svtx(ip a (x)-'&{x)) , (17) 

where we have written A(x) = \A(x)\e % ^ x \ Because 
|$+(0)| = \$-(L)\ = 1 for (w - Voo) 2 < A^, the initial 
values V?+(0) and ifi-(L) are real. Hence the solutions of 
Eq.(|l7|) remain real, which implies that |4> Q (x)| = 1 for 
all x. From Eq.(fLq) we obtain for the initial values 



tan^j + (0) = tan<yS_(L) = 



W — Va, 



-(lu- V^f 



(18) 



Note that the tp a (x) are unreduced phases which are not 
limited to take values between and 27r. In terms of the 
<f a (x) the local DOS can be written as 



p(x,u) 



.llmcotf^-^+iO 

7T V 2 



2E% = _ 00 6(<p + (x)-i P -(x)-2nm) . (19) 



We now make the w-dependence of ip a (x) explicit again. 
Since the right hand side of Eq. (O) is a 27r-periodic func- 
tion of (p a (x) it follows that if ip+Jx, oj) — ip-(x, ui) = 2nm 
for one x, this must also be true for all x. This implies 
that only for discrete values of lu does Eq.([l9]) yield a 
contribution to the local DOS. We get a delta-peak con- 
tribution to the total DOS if tp+(L, to) — 2irm + Lp + (Q, ui), 
where <p+(0, ui) = tp-(L,uj) is given by Eq.(|l8j). Since 
d ul tp + {x,u}) > 0, the integrated total DOS is given by 



1 

I 



<P+{L,u)) - <p+(L,0) 
2n 



C(u) 



(20) 



inl 
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where [z]i nt gives the integer value of z, and C(v) = 
(ip+(0,uj) — (fi+(0, 0))/2tt is a finite size correction of or- 
der unity that depends on the initial condition. For real 
A(x), a similar equation has been derived by Lifshits, 
Gredeskul, and Pastur m within the phase formalism. 
While these authors use a rather unphysical boundary 
condition, we can cope with arbitrary Aqo and Voo- In 
the thermodynamic limit the integrated DOS is indepen- 
dently of the boundary conditions given by 



A = A s g , A„+i = a„A„ + y/l - of A s g n +i , (24) 



AF(w) = lim (<p + (L,w) -<p+{L,0))/2irL 



(21) 



For large frequencies we recover the classical high- 
frequency limit Afa(uj) = w/tc, so that the DOS p(w) — 
dujAf(uj) is given by po — The white noise limit is 

also easily recovered: in this case Eq.(|l7|) implies that 
the probability distribution of <p+(x) satisfies a Fokker- 
Planck equation, which was first solved by Ovchinnikov 
and Erikhman [El for the commensurate case. For the 
most general case with complex A(x) see Ref. R|. 

Numerical algorithm - In the following we present an 
exact algorithm which allows to numerically calculate the 
(integrated) DOS for stepwise constant potentials. By 
choosing the step size sufficiently small, arbitrarily given 
potentials may be approximated in this way. Assuming 
that in the open intervals ]x n , x n +i[ the potentials A(x) 
and V(x) are given by the constants A„ and V n , the ma- 
trix U+(x) can be written as a finite product of matrices 
of the form 



-iM n 8 7l _ 



i " if; 



-i sinh[\/|A„| 2 - Q%d n 



= cosh[V|AnP ~ vl6 n ]a 
A n a + - A*cr_ + uj n a 3 



(22) 



where w„ = u) — V n and S n = Xjj+i — x n . For the Riccati 
variable & + (x) satisfying Eq.(Qq) this implies the recur- 
rence relation 



*+(Zn+l) 



r* n + t ra $ + (x„) 

t* n +T n ^+{x n ) 



(23) 



We found that for a given realization of A„ and V n it 
is easier to calculate the dynamics of $+(x) and to keep 
track of its phase than to directly solve Eq.([l7|). When- 
ever Re$+(x) > and there is a sign change in Im<I> + (x), 
the winding number [ip + (x)/2ir]- mt is changed by one. To 
detect all such changes we demand that the length S„ of 
all intervals satisfies 2(|u) n | + |A n |)5„ < ir/2. Since very 
long chains show a self- averaging effect, we only need to 
simulate one typical chain to obtain the average DOS. 

To generate Gaussian disorder with correlation length 
£ we have found the following realization of an Ornstein- 
Uhlenbeck process |i~9|, which is much simpler than the 
algorithm proposed in Ref. ||. Using the Box-Muller al- 
gorithm p0| , we generate independent Gaussian random 
numbers g n with (g n ) — and (g^) = 1. For real A(x) 
we set A n = A av + A„ and generate the A„ recursively 
according to 



where 



It is straightforward to show that 



this Markov process indeed leads to a Gaussian random 
process with the desired properties. Obvious advantages 
of our algorithm are that the random variables A n can 
be generated simultaneously with the iteration of the re- 
currence relation (p3), and that A n +i depends only on 
the previous A n , so that the implementation of this algo- 
rithm requires practically no memory space. Of course, 
our algorithm can also be used to generate V n , and in 
the complex case ReA„ and ImA„ can be generated by 
replacing A s by A s /v2. 

Results - In Fig.J^ we show our numerical results for 
p(u)/po for V(x) = A av = and real A(x). Except 
for A s £ = 1000, 0.2 we have chosen the same values 
of the dimensionless parameter A s £ as in Fig. 7 of Ref. 
fiof . One clearly sees the Dyson singularity, which ex- 
ists for any finite value of £ and overshadows the pseu- 
dogap at sufficiently small energies. This Dyson singu- 
larity is missed by Sadovskii's algorithm |L0]]. On the 
other hand, for complex A(x) this algorithm turns out 
to be qualitatively correct which can be seen by com- 
paring our data in Fig.^ with those in Fig. 5 of Ref. 
p0| . For a more quantitative comparison, the trian- 
gles (real A(x)) and diamonds (complex A(x)) in Fig.|| 
show the DOS p(uJ*) at the energy cu* where p(ui) as- 
sumes its minimum. Note that in the incommensurate 
case u>* — 0. The numerical errors (which are mainly 
due to the finite length of the chain) are smaller than 
the size of the symbols. The straight lines are fits to 
power-laws p(uj*)/po = A(A S £) _A1 . For real A(x) we 
obtain A = 0.482 ± 0.010, p = 0.3526 ± 0.0043, while 
for complex A(x) our result is A — 0.6397 ± 0.0066 and 
p = 0.6397 ± 0.0024, i.e. within numerical accuracy we 
find A = p. The circles in Fig.|| show for real A(x) the 
energy scale ui* where p{ui) is minimal. The long solid 
line is a fit to a power-law lu*/A s — £?(A S £) -7 , with 
B = 0.2931 ± 0.0074 and 7 = 0.3513 ± 0.0051. Here 
we find within numerical accuracy p = 7. The propor- 
tionality of p(co*) to the energy scale u>*, which can be 
interpreted as the width of the Dyson singularity, can 
also directly be seen in Fig.0. Finally we note that for 
A s £ ;$ 0.2 our algorithm produces results consistent with 
the white noise limit A s £ <C 1. While p(0) — > 1 in the in- 
commensurate case, in the commensurate case we obtain 
from the exact solution of Ovchinnikov and Erikhman j^] 
p{uj*)/p -> 0.9636, and uj* -> 1.2514Af£, which deter- 
mines the short solid line in Fig.|| describing cj*(£) in the 
white-noise limit. 

Summary - We have developed a powerful numerical 
algorithm to calculate the average DOS of the FGM with 
very high accuracy. The algorithm can be used for arbi- 
trary forward and backward scattering potentials, so that 
it is not restricted to the case of vanishing averages and 
Gaussian statistics which we further considered in this 
work. Our main results are: (a) for commensurate chains 
in the absence of forward scattering the DOS exhibits for 
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large £ a pseudogap and a Dyson singularity. We have 
explicitly calculated the width of the Dyson singularity 
as a function of £. The most promising experimental sys- 
tems to detect Dyson singularities are spin chains Ijj . (b) 
In the incommensurate case the algorithm proposed by 
Sadovskii |L(| is qualitatively correct. However, his result 
(p(0)) oc ^ 1 i 2 is incorrect. This should be kept in mind 
for a quantitative comparison between experimental data 
PI and calculations based on Sadovskii's algorithm. 

We thank K. Schonhammer for discussions. This work 
was financially supported by the DFG (Grants No. Ko 
1442/3-1 and Ko 1442/4-1). 
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FIG. 1. Average DOS for real A(x) with A a L = 10 7 , 
V(x) = A av = 0, and A s £ = 1000, 100, 10, 2, 1, 0.5, 0.2. 
The minimal DOS p(co*) decreases with increasing A s £. 




FIG. 2. Average DOS for complex A(x) with A S L = 10 7 , 
V(x) = A av = 0, and A s £ = 1000, 100, 10, 2, 1, 0.5, 0.2. p(0) 
decreases with increasing A s £. 



1.00 



o 0.10 : 



^3, 



0.01 




1.00 10.00 



FIG. 3. Double-logarithmic plot of p(u>*)/po as a function 
of l/Asf for real A(x) (triangles) and complex A(x) (dia- 
monds), where u>* is the energy for which the DOS assumes 
its minimum. While uj* = for complex A(x), the circles 
give the double-logarithmic plot of ui* / A s for real A(x) as a 
function of 1/A S £. 
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